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PARAMETER-ROBUST DISCRETIZATION AND 
PRECONDITIONING OF BIOT’S CONSOLIDATION MODEL 

JEONGHUN J. LEE, KENT-ANDRE MARDAL AND RAGNAR WINTHER 


Abstract. Biot’s consolidation model in poroelasticity has a number of appli¬ 
cations in science, medicine, and engineering. The model depends on various 
parameters, and in practical applications these parameters ranges over several 
orders of magnitude. A current challenge is to design discretization techniques 
and solution algorithms that are well behaved with respect to these varia¬ 
tions. The purpose of this paper is to study finite element discretizations of 
this model and construct block diagonal preconditioners for the discrete Biot 
systems. The approach taken here is to consider the stability of the problem 
in non-standard or weighted Hilbert spaces and employ the operator precon¬ 
ditioning approach. We derive preconditioners that are robust with respect to 
both the variations of the parameters and the mesh refinement. The param¬ 
eters of interest are small time-step sizes, large bulk and shear moduli, and 
small hydraulic conductivity. 


1. Introduction 

Biot’s consolidation model describes the deformation of an elastic porous medium 
and the viscous fluid flow inside when the porous medium is saturated by the fluid. 
The unknowns are the displacement of the elastic medium, u 1 and the fluid pres¬ 
sure, pf • In homogeneous isotropic linear elastic porous media, the equations for 
the quasi-static Biot model are: 

— div(2 ge(u) + A div ul_ — apFl.) = f, 
sqPf + ol div ii — di v(kVpf) = g , 

where the dots represent time derivatives, g and A are the Lame coefficients of 
clastic medium, e(u ) is the symmetric gradient of u, J is the n x n identity matrix, 
So > 0 is the constrained specific storage coefficient, n > 0 is the hydraulic con¬ 
ductivity determined by the permeability of medium and the fluid viscosity, and 
a > 0 is the Biot-Willis constant which is close to 1. The system m can be posed 
on a bounded domain in two and three space dimensions. The given functions f 
and g represent body force and source/sink of fluid, respectively. We will assume 
throughout the paper that the parameters g, A, and So are scalar functions on the 
domain, while in general k can be a symmetric positive definite matrix-valued func¬ 
tion. To be a well-posed mathematical problem, the system OD needs appropriate 
boundary and initial conditions. A discussion of general boundary conditions for 
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the Biot system will be given in Remark 13.31 below. Furthermore, a mathematical 
discussion of well-posedness of this model can be found in |T|. 

Due to importance of Biot’s model in applications, ranging from geoscience 
to medicine, finite element methods for the model have been studied by many 
researchers. For example, various primal methods are studied in mixed 

methods in 00 Cl, Galerkin least square methods in [8], discontinuous Galerkin 
methods in [9j, and combinations of different methods in pzanumm, but this 
list is by no means complete. 

It is important to construct numerical methods which are robust with respect 
to variation of model parameters since this variation in many practical problems 
is quite large. For example, relevant parameters in the soft tissue of the central 
nervous system are Young’s modulus of 1 — 60 kPa, Poisson ratio from 0.3 to 
almost 0.5 (0.499 in [15), and the permeability is 10 -14 — 10 _16 m 2 [HI 116... In 
geophysics, Young’s modulus is typically in the order of GPa, Poisson ratio 0.1 —0.3, 
while the permeability may vary from approximately 10~ 9 to 10 _21 ?ti 2 [miii]. 
Relations of Young’s modulus E, Poisson ratio v and the two elastic moduli /r, 
A are n = E/2(l + v) and A = Ev/{\ + 1 — 2v). Consequently, /i and A are 

in the ranges of 300 — 500 MPa and 100 — 500 MPa, respectively, in geoscience 
applications, whereas corresponding numbers are /r and A in the ranges 300 — 2000 
Pa and 500 — 10 6 Pa in neurological applications. 

However, in the present paper we will not limit ourselves to the study of ro¬ 
bustness with respect to model parameters of the finite element discretization of 
Biot’s model. In fact, our main concern is to be able to construct preconditioners 
for the discrete systems which are well behaved both with respect to variations of 
the model parameters and the refinement of the discretization. When large discrete 
systems are solved by iterative methods, the convergence rate depends heavily on 
the construction of suitable preconditioners. Such preconditioners for finite element 
discretizations of Biot’s model have been studied by many authors, cf. for example 
[m EnumEn [ 23 ] . Recently, there is also an emerging interest for preconditioners 
which are robust with respect to model parameters [23 121] • However, robustness 
with respect to all model parameters remains challenging. In particular, we will 
derive preconditioners that are robust as the medium approaches the incompress¬ 
ibility limit while the permeability is low. In our experience this represents the most 
difficult case, and it is also the case that occurs in many biomechanical applications. 

The purpose of this paper is to develop a stable finite element method for Biot’s 
model, and a corresponding preconditioner for the associated discrete systems, such 
that the preconditioned systems have condition numbers which are robust with re¬ 
spect to variations of model parameters. More precisely, we aim to have a precon¬ 
ditioned system which is robust for small k, small time-steps, large A, large /i, and 
mesh refinements. In order to obtain such a parameter-robust preconditioner we 
employ the operator preconditioning framework of [55] • It turns out that typical 
formulations of Biot’s model are not appropriate to apply the framework, so we de¬ 
velop a new three-field formulation of Biot’s model and propose a parameter-robust 
block diagonal preconditioner for it. 

The present paper is organized as follows. In Section 2, we introduce some 
notation and conventions that will be used throughout the paper. Furthermore, 
we briefly discuss the preconditioning framework of [25] based on parameter-robust 
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stability of the continuous problems, and illustrate this with some numerical exam¬ 
ples based on simplified models which can be seen as subsystems of the Biot system. 
In Section 3 we explain some difficulties related to more common formulations of 
the Biot system, and as a consequence we motivate a new three-field formulation. 
The discussion of finite element discretizations based on this three-field formulation 
is given in Section 4, and the stability results are used to motivate the construction 
of parameter robust preconditioners. The implementation of a special operator re¬ 
lated to one of the blocks a block diagonal preconditioner is discussed in Section 5. 
Finally, in Section 6 we present some numerical experiments which illustrate our 
theoretical results. 


2. Preliminaries 

The system HU can in principle be studied on rather general domains Q in two 
or three dimensions. However, our main goal is to study finite element approxima¬ 
tions of this system, and therefore we will assume throughout this paper that Q is 
a bounded polyhedral domain in R ra , with n = 2 or 3. We will use H k = H k (Q) 
to denote the Sobolev space of functions on H with k derivatives in L 2 and the 
corresponding norm is denoted by || • ||a-. Further, let Hq be the closure of Co°(fl) 
in H k with dual space denoted by H~ k and (•,•) denote the L 2 inner product of 
scalar, vector, and matrix valued functions as well as the duality paring between Hq 
and H~ k . The space L\ is the space of L 2 functions with mean value zero. Bold¬ 
face symbols are used to denote vector valued functions or spaces, and symbols of 
boldface with underline are used to denote matrix valued functions. 

Throughout this paper we use A < B to denote the inequality A < CB with a 
generic constant C > 0 which is independent of the discretization parameters and 
the model parameters, and A ~ B will stand for A < B and B < A. If needed, we 
will use C to denote generic positive constants in inequalities. For a scalar valued 
function g , Vg is a (column) vector valued function. For a matrix valued function 
g , div is understood as a row-wise divergence which results in the vector valued 
function div < 7 . Adopting these conventions, the equations (ED are well-defined. 

2.1. Preconditioning of parameter-dependent systems. Let us briefly review 
the abstract framework of parameter-robust preconditioning in [23] ■ Let X be a 
separable, real Hilbert space with inner product {■r)x and the associated norm 
|| • ||x- For two Hilbert spaces X and Y. C(X,Y) is the Hilbert space of bounded 
linear maps from X to Y. Let us denote the dual space of X by X* and the duality 
pairing of X and X* by (■,•). Suppose that A £ C(X,X*) is invertible and also 
symmetric in the sense that 

(Ax,y) = (x,Ay), x,y € X. 

For given / £ X* consider a problem finding x £ X such that 
(2.1) Ax = f. 

Its preconditioned problem with a symmetric isomorphism B £ C{X*,X) is 

BAx = Bf. 

The convergence rate of a Krylov space method for this problem can be bounded 
by the condition number, K(BA ), given by 

K(BA) := \\BA\\c(x,x)\\(BA)- 1 \\ C (x,xy 


4 


JEONGHUN J. LEE, KENT-ANDRE MARDAL AND RAGNAR WINTHER 


Parameter-dependent problems are handled in this framework as follows. Let e 
denote a collection of parameters, and A e the parameter-dependent coefficient op¬ 
erator. A systematic way to construct an e-robust preconditioner B e , as proposed in 
|25| , is to consider the mapping property of A e in e-dependent Hilbert spaces X e and 
X*. The key property is to choose the spaces X e and X* such that A e is a map from 
X e to X*, and with corresponding operator norms ||»4 e ||£(x e ,X|) and || Aj 1 \\c(x*,x e ) 
bounded independently of e. In this case the preferred preconditioner, B e , is a 
map from X* to X e with the property that ||# e ||,c(x*,x s ) an d ||B“ 1 ||£(x e ,x») are 
bounded independently of e. If such an operator B e is identified, then the condition 
number K(B e A e ) is bounded independently of e, since both B £ A E and {B e A e )~ l 
are operators on C(X e , X e ), with corresponding operator norms bounded indepen¬ 
dently of e. 

As we will illustrate below the discussion outlined above can often most easily 
be done in the continuous setting. On the other hand, in a computational setting 
we need preconditioners for the corresponding discrete problems. In fact, if we 
utilize a finite element discretization which is uniformly stable with respect to the 
parameters, then the structure of the preconditioners in the discrete case will be the 
natural discrete analogs of the preconditioners in the continuous case. However, the 
preconditioners derived by the procedure above will often require exact inverses of 
operators which cannot be inverted cheaply. Therefore, in order to obtain effective 
robust preconditioners in the discrete case, we also have to replace these exact 
inverses by related equivalent operators, often obtained by common procedures 
such as multilevel methods or domain decomposition methods. We refer to [25] 
and the examples below for more details. 

A challenge of the Biot system is the dependency of several different and inde¬ 
pendent parameters like the Lame elastic parameters as well as parameters related 
to porous flow such as permeability and the Biot-Willis constant. The aim is to 
achieve robustness with respect to all model parameters, as well as the resolution of 
the discretization. To motivate this discussion we start by considering two simpli¬ 
fied examples related to the Biot equations. The first example illustrates the case 
where the permeability tends to zero, while the Lame parameters are of unit scale. 
In the second example we consider the case where the elastic material tends towards 
the incompressible limit. Both examples are special cases of the Biot system. 

Example 2.1. Consider a system of equations 

—Am + Vp = /, 

( 2 . 2 ) 

— divM + div(reVp) = g, 

with unknowns u : H — > R” and p : H — > R. As boundary conditions we use 
homogeneous Dirichlet condition for m, i.e., u\ga = 0, while we use homogeneous 
Neumann condition for p. The parameter re > 0 is taken to be a constant in this 
example. A variational formulation of this problem is to find (u,p) £ Hi x H 1 flLg 


satisfying 


(Vm, Vv) - (p,divv) = (f,v), Vv £ Hi, 

— (div m, q) - (reVp,Vg) = {g,q), Wq £ H 1 n L 2 0 

This system has a form of (12.11) with X = Hi x H 1 D Lq an d 
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If we define X K as the Hilbert space H q x H 1 n L\ with ^-dependent norm given 
by 


ll(«»P)llx„ = Hli + Ibllo + «l|Vp||g, 


then one can check that A : X K —> X* is an isomorphism with corresponding 
operator norms of A and A~ l bounded independently of n. Here the norm on 
X* D L 2 x L 2 is defined from the norm on X K by extending the L 2 inner product 
to a duality pairing. To define a robust preconditioner we need to identify an 
isomorphism B : X* —> X K with corresponding operator norms of B and B~ x 
bounded independently of k. A natural choice is a block-diagonal operator of the 
form 


(2.3) 


B = 


-A- 1 0 

0 (/ — kA) _1 


However, the operator B only indicates the structure of the desired preconditioner of 
the discrete system. We will discretize this problem on the unit square in R 2 by the 
lowest order Taylor-Hood element with respect to a mesh of uniform squares. This 
method is uniformly stable with respect to n in the proper norms. Furthermore, 
to obtain an effective preconditioner the exact inverses in the definition of B are 
replaced by corresponding algebraic multigrid preconditioners for the operators — A 
and I — kA. implemented in the software library Hypre [2H] with default settings. 
The preconditioned system is implemented using cbc.block [27] and FEniCS [251 . 
The eigenvalue estimates are obtained by the conjugate gradient method of normal 
equation of the system with convergence criterion 10~ , and we refer to (27l for 
more details. The same setup is used throughout the paper. 

We present numerical results in TablcQ] The numbers of iterations and condition 
numbers increase as the value of k decreases but they are asymptotically stable and 
are still bounded in the limit case n = 0, which is the Stokes equation. We remark 
that the zero eigenvalue of the system associated with H 1 D Lq is ignored in the 
computation of condition numbers. 


TABLE 1. Number of iterations of MinRes solver of system (12.21) with 
algebraic multigrid (AMG) preconditioner of the structure (12.31) . Esti¬ 
mates of the condition numbers of the preconditioned system are given 
in parenthesis. (H = unit square, partitioned as bisections of N x N 
rectangles, convergence criteriorQ with relative residual of 10 -6 ) 



16 

32 

N 

64 

128 

256 

10° 

13 (1.3) 

13 (1.2) 

14 (1.2) 

14 (1.2) 

14 (1.2) 

nr 1 

16 (1.7) 

16 (1.8) 

16 (1.8) 

16 (1.8) 

16 (1.8) 

icr 2 

22 (3.1) 

24 (3.4) 

23 (3.6) 

23 (3.8) 

24 (3.9) 

lo -3 

30 (4.6) 

29 (4.9) 

29 (5.4) 

30 (5.8) 

30 (6.1) 

K io- 4 

35 (6.1) 

36 (6.2) 

36 (6.4) 

35 (6.6) 

34 (7.0) 

icr 5 

38 (7.0) 

39 (7.8) 

41 (7.9) 

39 (7.7) 

39 (7.7) 

icr 6 

38 (7.1) 

40 (8.3) 

42 (9.2) 

44 (9.4) 

43 (9.0) 

0 

38 (7.1) 

42 (8.4) 

44 (9.5) 

47 (10.5) 

48 (11.2) 
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Example 2.2. Consider a system related to the Lame problem in linear elasticity: 
Find u : fl -A R”, p : —> M for 

- dive (it) - Vp = /, 

( 2 - 4 ) 1 

dlV U ~ J P = 9 ’ 

with u\qq = 0, where 1 < A < +oo is a positive constant and e(u) is the symmetric 
gradient of u. When g = 0, this is the Lame problem and p = A div u is called 
“solid pressure”. Its variational form is to find u G Hq and p £ L 2 such that 

(e(u),e(v)) + (p, div v) = (f,v), Vn G Hi, 

(di .vu,q) - ~{p,q) =g, \/qeL 2 . 

This is a saddle point problem with a stabilizing term — (1/A)(p, q), and the sta¬ 
bilization becomes weaker as A becomes larger. Let X = Hi x L 2 be the Hilbert 
space with standard norm. In the limit when A = +oo the system is not stable 
in these norms. This is due to the fact that div Hi C L 2 and as a consequence 
the Brezzi condition for stability of saddle point problem is not fulfilled [29]. More 
precisely, div Hi can control only the L 2 norm of the mean-value zero part of p , 
and the stabilizing term is needed to control the mean-value part of p. Since the 
stabilizing term is dependent on A, we need a A-dependent norm on Hi x L 2 to 
have A-independent stability of the system. 

Before we define an appropriate A-dependent norm, we need some preliminaries. 
Let P m be the linear operator in L 2 such that 

Pm< ^ := {w j n ^ dx ) Xn ’ Vcfr&L 2 , 

where yo is the characteristic function on H and |f2| is the Lebesgue measure of fl. 
Notice that P m </> and <fi— P m 4> are the decomposition of <j) into its mean-value part 
and mean-value zero part. For q G L 2 we denote its mean-value and mean-value 
zero parts by 

(2.5) qm := Pmq, qo = q qm- 

We now define a Hilbert space X\ by 

IIK9)IIa a = Mi + Qlkm||o + IMIo) » 


then the system is A-independent stable in X\. We will not give a detailed proof of 
stability here since it can be obtained by modifying the proof of Theorem 13.21 below. 
But the rough explanation is that the mean value of p cannot be control by the 
inf-sup condition, and therefore this part of the norm has to be weighted properly 
in balance with the stabilizing term A -1 ||p m ||g, while the rest of p is controlled by 
the inf-sup condition. 


We present numerical results for two different preconditioners, 
space X and X\ lead to preconditioners of the forms 


The Hilbert 


( 2 . 6 ) 


Pi = 


-A- 

0 


0 

I- 1 


B 2 = 


-A- 

0 


(lo 


0 

1 / 

\ ±rT 


■^Convergence criterion is (Br ^, rp/(Br^. tq) < 10 6 where B is preconditioner and rf ; . is the 
residual of k- th iteration. 
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Here the appearance of the operator 7 _1 calls for an explanation. In fact, the 
operator / should not be thought of as the identity operator on the Hilbert space 
L 2 , but rather as the Riesz map between this space and its dual. In particular, in 
the corresponding discrete setting the operator / -1 is typically represented by the 
inverse of a mass matrix. In a similar manner, the operators I o and I m should be 
interpreted as maps of L 2 into the duals of L\ and its complement. We refer to [25] 
Section 6] for more details. 

We employ the lowest order Taylor-Hood discretization and test the efficiency of 
the preconditioners on different refinements of the unit square. The preconditioner 
Hi is implemented by replacing the exact inverse of — A by an algebraic multigrid 
operator, and we use the Jacobi iteration to approximate / . The construction 

of H 2 is technical due to the second block, and we postpone the details of the 
construction to Section [5] 

Numerical results for the preconditioners B\ and B 2 are given in Tables [2] and 
[3] The results for the preconditioner with structure B 2 is completely satisfactory. 
Both the number of iterations and the condition numbers appear to be uniformly 
bounded with respect to A and mesh refinement. However, the results for the 
preconditioner B\ are different. In this case the condition numbers appear to grow 
linearly with A, while the number of iterations still seems to be bounded, even if they 
appear to be slightly less robust in this case. So in the case of the preconditioner B\ 
the condition number will not lead to a sharp bound on the number of iterations. 
In fact, by comparing the operator B\ with the uniform preconditioner B 2 one 
can argue that the operator B\A has one isolated eigenvalue that tends to 0 as A 
increases, while the rest of the spectrum lie on an interval bounded independently 
of the mesh resolution and A. In cases where only few eigenvalues are outside a 
bounded spectrum, it is well-known that the Conjugate Gradient and the Minimum 
Residual methods are very efficient [30] J31] and therefore B\ is about as efficient as 
B 2 , but slightly less robust. 


TABLE 2. Number of iterations of MinRes solver of system (12.41) with 
preconditioner of the form B\ in (12.61) . Estimates of the condition num¬ 
bers of the preconditioned system are given in parenthesis, (fi = unit 
square, partitioned as bisections of N x N rectangles, convergence cri¬ 
terion with relative residual of 1CP 6 ) 



16 

32 

N 

64 

128 

256 

10° 

29 (3.5) 

29 (3.6) 

29 (3.6) 

29 (3.6) 

29 (3.6) 

10 1 

40 (10.8) 

41 (10.8) 

38 (10.9) 

38 (10.9) 

36 (10.9) 

10 2 

53 (95.6) 

59 (96.2) 

54 (96.7) 

53 (96.9) 

52 (97.0) 

A 10 3 

60 (958) 

62 (964) 

62 (969) 

61 (972) 

60 (973) 

10 4 

66 (9589) 

69 (9649) 

44 (9697) 

43 (9724) 

42 (9734) 

10 5 

44 (95892) 

44 (96501) 

44 (96981) 

43 (97248) 

42 (97344) 

10 6 

44 (958942) 

44 (964970) 

44 (969940) 

43 (972566) 

41 (973587) 


2.2. Parameter rescaling of the Biot system. The systems in the above two 
examples have only one parameter, so it is relatively easy to find function spaces 
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TABLE 3. Number of iterations of MinRes solver of system (12.41) with 
preconditioner of the form B 2 in (12.61) . Estimates of the condition num¬ 
bers of the preconditioned system are given in parenthesis, (Cl = unit 
square, partitioned as bisections of N x N rectangles, convergence cri¬ 
terion with relative residual of 1CP 6 ) 



16 

32 

N 

64 

128 

256 

10° 

29 (3.5) 

29 (3.6) 

29 (3.6) 

29 (3.6) 

29 (3.6) 

10 1 

41 (11.6) 

41 (11.6) 

41 (11.6) 

38 (11.6) 

38 (11.6) 

10 2 

46 (18.4) 

45 (18.4) 

44 (18.4) 

44 (18.4) 

44 (18.4) 

A 10 3 

46 (19.6) 

46 (19.6) 

45 (19.6) 

45 (19.6) 

44 (19.6) 

10 4 

45 (19.7) 

46 (19.7) 

45 (19.7) 

45 (19.7) 

44 (19.7) 

10 5 

45 (19.7) 

46 (19.7) 

46 (19.7) 

44 (19.7) 

44 (19.7) 

10 6 

46 (19.7) 

46 (19.7) 

43 (19.7) 

44 (19.7) 

44 (19.7) 


and parameter-dependent norms such that the aforementioned preconditioner con¬ 
struction is applicable. However, the Biot system has several parameters of different 
ranges, so it is not easy to find appropriate function spaces and their parameter- 
dependent norms. In the rest of this section, we will rescale parameters of the Biot 
system and reduce it to a problem with three parameters. This procedure will not 
only simplify the problem but also clarify intrinsic parameters of the system. We 
emphasize that p, A, so, and n are allowed to be functions on the domain, while a 
is assumed to be a constant. 

Recall that when we solve a time-dependent problem numerically, we discretize 
the problem in time and solve a static problem at each time step, so preconditioning 
of time-dependent problem is reduced to preconditioning of static problem at each 
time step. If we consider an implicit time discretization (e.g., the backward Euler 
method) applied to (11.11) with time-step size 5 2 (0 < <5 < 1), and multiply the 
second equation with — 5 2 , then we obtain a static problem 

— div(2 pe(u) + A div wJ — cxpfL ) = /, 

—sqPf — ot div u + 5 2 di v(k\7pf) = g 
with some right hand side g. 

To reduce this system further, we recall the physical derivation of so- The storage 
coefficient so is the increase of the amount of fluid for the unit increase of fluid 
pressure, when volumetric strain is kept as constant. More precisely, so = 4>cf + 
(1 — 4>)cs where <j> is the porosity of solid, and cf > 0, cs > 0 are compressibilities 
of the fluid and solid. For derivation of these equations from physical modeling, 
we refer to standard porous media references, for instance, [32] . The parameter cs 
and other parameters p and A are related so that cs ~ 1/(2 p/n + A) holds with 
n, the spatial dimension of Cl. In many practical applications, p < A and cf ~ 0 
hold, so we have So ~ 1/A. Furthermore, a is a constant close to 1, so we will 
assume that So scales like a 2 /X. Therefore, to reduce the number of parameters in 
our system we will simply let s 0 = a 2 /A for the rest of the discussion in this paper. 
We emphasize that the equality is not essential. Our analysis below can easily be 
adopted to the situation where so scales like o? /A. However, this rather artificial 
expression of sq is useful when we normalize the system later. In addition, S 2 k can 
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be regarded as one small parameter because the hydraulic conductivity k is small 
in general. As a consequence, introducing k = S 2 k , we have a simplified system 


— div(2 pe(u) + A div wJ — cxpfL ) = /, 
a 2 

— ~pf — ol div u + di v(k\7pf) = g, 

A 

with 1 < A < +oo, small k, and a ~ 1. However, in practical applications, p can be 
much larger than a, so we rescale the above equations to include this factor. To do 
so we assume p is a spatial function with a uniform scale, i.e., there is a constant 
p such that p/p ~ 1. Let 


A' = A 


/ M 

g = -, 


a = 


a 

2/2 ’ r p ' 2 / 2 ’ 

Dividing the above two equations by 2/2, the final simplified equations are 


, _ « 

2 / 2 ’ ^ 2 / 2 ’ 


_ a 

9 2 / 2 ' 


(2.7) 


— di v(p'e(u) — A' div «/ — a'ppL) = f - 


(qQ 2 

A' 


p F - 


a' div it + div(K , Vpj?) = g', 


where a' is a small positive constant, X' and p' are positive scalar functions such 
that X' is bounded from below and p' ~ 1, while n' is a positive definite matrix 
valued function with eigenvalues bounded from above. 

In summary, we have reduced the original Biot system to a system with three 
intrinsic parameters A', a', and k' which all may be unbounded. More precisely, 
the positive constant a 1 may be arbitrarily small, the scalar function X' may be 
arbitrarily large, while the positive eigenvalues of k' may be close to zero. Since 
p’ is bounded from above and below this parameter has no essential effect on the 
properties of the system. Therefore, to reduce the number of parameters of the 
system, we take p’ = 1 in the discussion below. 

For the biomedical applications discussed in the introduction in units of Pascal, 
gram, milimeter and second, p was in the order of 1 — 60 kPa which makes X', a ', «/ 
in the ranges of 0.25 — 500, 10 -3 — 10 -5 , 10 -8 — 10 -12 , respectively. In geoscience, 
a representative p is 10 GPa and units for pressure, viscosity and permeability are 
pounds per square inch (psi), centi Poise (cP), mili Darcy (rnD). Using these units, 
A', a', k' are in the ranges of 0.25 — 3.5, 10 -10 , 10~ 9 — 10 -13 , respectively. 


3. Parameter-robust stability of the continuous problems 

For the rest of this paper we will discuss a system of the form (12.71) . where the 
parameter p! = 1. By omitting the primes on the parameters we obtain a system 
of the form 

— div(e(it) — A div ul_ — upfL) = f , 

(3-1) a 2 

— ~pf — ol div u + di v(kX7pf) = g■ 

A 

Throughout the rest of the paper we will assume that the parameters A, a, and k 
satisfy 

(3.2) 1 < A < + 00 , 0 < a < 1, 0 < k < 1, 

where the assumption on the matrix valued function k has the interpretation that 
k is uniformly positive definite, and with all eigenvalues bounded above by one. 
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3.1. Difficulties in typical formulations. In the discussion of the two examples 
in Section 2 we saw that simplified versions of the Biot system were efficiently han¬ 
dled with straightforward extensions of standard preconditioners for Stokes prob¬ 
lem. In particular, in Examnlc l2.il the presence of a small permeability was handled 
by extending a standard Stokes preconditioner in the canonical way with an oper¬ 
ator of the form V ■ («V) applied to the pressure. Furthermore, in Example 12.21 
the use of “solid pressure” gave a stable formulation and an efficient preconditioner 
even in the limit of an incompressible material. We will now demonstrate that 
extending these two approaches to the Biot system is not straightforward. 

For simplicity we will first consider the system with homogeneous Dirichlet 
boundary conditions of the form u = 0 and p F = 0 on dfl. A discussion of 
other possible boundary conditions is given in Remark 13.31 below. Throughout this 
subsection we will make the simplifying assumption that A is a constant, and we 
will illustrate that parameter-robust preconditioning is difficult even in that case. 
For a variational formulation of (13.11) with Dirichlet boundary conditions, we will 
use the function spaces Hq and FFq for the unknowns u and p F , respectively, and 
obtain 


(e(u), e(i>)) + A(divit, divu) — a(p F , divu) = 
a 2 

—a(div u 1 q F ) - —( p Fl q F ) - ( K.X7p F ,X7q F ) = {g,q F ), 


Vw G Hi, 
\/q F G Hq . 


In matrix form, the system is 


(3.3) 



/— div(e + A div J) 
\ —a div 



aV 

+ div(h'.V) 




This system has a perturbed saddle point problem structure. Following the precon¬ 
ditioner construction framework in the previous section, we define Hilbert spaces 
Vi and Qp,i with parameter-dependent norms on FTq and Hq, 

IMIvg : = kMllo + A ll divu||o, v G Hi, 


\\Qf\\q f<1 ■= — IkiHIo + {n\/q F ,Vq F ), 


q F G Flo 1 . 


Then we are able to show that A : V\ x Qf,i —> V\ x Q F1 in (13.31) is an iso¬ 
morphism with upper and lower bounds uniform in A, k, and a. However, this 
formulation has a nontrivial difficulty to achieve parameter-robust preconditioner 
for its discrete counterpart. For example, if we consider a block-diagonal precon¬ 
ditioner as in examples in the previous section, we need a good preconditioner of 
— dive — A grad div in first block of the preconditioner. However, usually accepted 
preconditioners (e.g., algebraic multigrid preconditioner) do not perform well when 
A is large. This is observed in [23] for the simplified McKenzie equations, and 
can be explained by the fact that it is hard to construct discretizations which are 
uniformly stable with respect to A. 

This difficulty is similar to volumetric locking problem in linear elasticity [ 33 ], 
which arises when A is very large. Thus, we expect that resolutions of the lock¬ 
ing problem in elasticity are useful to circumvent this preconditioning problem. 
There are two mathematically equivalent ways to avoid the locking problem: one 
is reduced integration technique [34] and the other is the mixed method (see, e.g., 
[551 l36l l37l 38] 1. However, both of them have technical difficulties in parameter- 
robust preconditioning. Here we will discuss the difficulty with the mixed approach. 
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Motivated by Example 12.21 and the mixed finite element technique to avoid the 
locking problem in linear elasticity, it is tempting to employ the “solid pressure” 
ps := —A div it. Recall that our purpose here is to show that this formulation is not 
appropriate for parameter-robust preconditioner construction, so we show the lack 
of stability and bad numerical results only for one specific case a = 1. Introducing 
ps = — A div it in the first equation of m, with some algebraic manipulation, we 
have a three-field formulation 

- div e(it) + Vps + VpF = f, 

(3.4) — div it — —ps = 0, 

A 

— div it — ~pf + div(KVp F ) = g■ 

A 

Since u £ H q, div it is mean-value zero, and therefore ps = —A div it is mean- 
value zero as well. This means that, Lq is an appropriate function space for ps in 
variational formulation. Thus, a variational form of ( l 3 ~ fl ) is to find ( u,ps,pf ) £ 
Hq x Lq x Hq such that 

(e(it),e(i;)) - (p s + p F , div v) = (f,v), Vv £ H\, 

— (div it, q s ) — ^iPSi Qs) = 0, Vq s £ L 2 0 , 

— (div it, q F ) - j(p F ,qF) - (kVpf, Vg F ) = (s,9f), £ H]. 

In order to have parameter-robust stability, we need to find parameter-dependent 
norms of Hq x Lq x Hq such that the above system gives a linear isomorphism 
from the Hilbert space to its dual space, and norms of the linear isomorphism and 
its inverse are independent of the parameters. The bilinear forms in the system, 


(3.5) 

(e(w),e( v)), 

u,v £ Hi, 

(3.6) 

(div v,p s ), 

v £ Hi, p s £ Lq, 

(3.7) 

(div v, p F ), 

v £ H l Q,p F £ Hi, 


j(ps,qs), 

Ps, Qs £ L 2 q, 

(3.8) 

j(pf,Qf) - («Vp F , Vg F ), 

Pf, Qf £ Hi, 


have to be bounded for the parameter-dependent norms, so some necessary condi¬ 
tions of the norms will be given. The bilinear form (13.51) suggests KI 1 -norm for JLj. 
To make (13.61) and (13.71) bounded, the chosen norms of Lq and Hq have to bound 
the standard L 2 -norms of ps and p F , respectively. Finally, (13.81) enforces the norm 
of Hq to be an upper bound of («: Vp F , Vpf) 1 / 2 . Thus the smallest possible norms 
for Hq, Lq, Hq from these observations are 

(3-9) ||e(u)|| 0 , Ibsllo, (||p F ||o + («;Vp F , Vp F )) 5 

for u £ Hi, p s £ L 2 q, p F £ Hi. 

We use V 2 , Qs, Qf ,2 to denote the Hilbert spaces on H J, Lq, H^ with the 
above norms. It is easy to check that all bilinear forms are bounded with these 
norms uniformly in A and k. In other word, the linear map from V 2 x Qs x Qf, 2 
to V 2 x Q*s x <5 F2 ; given by the above three-field formulation has a uniform bound 
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independent of the parameters. Unfortunately, it does not seem to be the case for 
the inverse of the linear map. Although the system is a stabilized saddle point 
problem, the stabilization terms —A ~ 1 (ps,Qs) and —A ~ 1 (p F ,q F ) are not enough 
to control the L 2 norms of ps and pf when A is very large. Therefore, we need to 
control the norms by v £ V 2 with inf-sup condition. In other words, we need 


(3.10) 


inf sup 

{qs,< if)£QsxQf ,2 D6V2 


(divt),g s +q F ) 

Hlidksllo + Mlo) 


>/?> 0 , 


with a constant f3 which is independent of the parameters. However, both qs and q F 
interact with div v in the bilinear form (div v, qs + q F ), and it is difficult to control 
two independent terms with one object, div v. When k is not small the stabilization 
term ( nVp Fl 'S/p F ) can be used to control the L 2 -norm of p F . However, as we have 
seen in the model reduction in the previous section, the smallness of n is given not 
only by small hydraulic conductivity in the model, but also by a small time-step. 
Thus it is inevitable to assume that k is small when we solve static problems at 
each time step. 


Example 3.1. We present a computational result which provides numerical evi¬ 
dence for the above discussion. Suppose that H is the unit square in R 2 and 

T d = {(x,y)eR 2 : (x, y) <E <9H, x < 1}. 

For simplicity of implementation, we assume that u is vanishing on U, not on 
dfl, and therefore the appropriate function space for ps is L 2 , since Adivu ^ Lq. 
This is a reasonable assumption since robust preconditioners should cover problems 
with general boundary conditions. For discretization we use elements inspired by 


TABLE 4. Number of iterations for different A and k with the precondi¬ 
tioner of the form in (13.1111 . (fl = unit square, partitioned as bisections 
of NxN rectangles, convergence criterion with relative residual of 1CF 6 ) 






k 

(k = 

10 fe ) 



N 

X 

0 

-1 

-2 

-3 

-4 

-5 

-6 


10° 

21 

20 

21 

22 

24 

25 

25 


10 2 

45 

46 

49 

88 

115 

111 

84 


10 4 

46 

46 

50 

102 

247 

359 

198 


10 6 

46 

46 

53 

102 

251 

451 

238 


10° 

20 

20 

20 

21 

23 

23 

24 


10 2 

44 

44 

49 

85 

107 

107 

92 


10 4 

44 

44 

48 

98 

242 

315 

207 


10 6 

44 

44 

50 

98 

244 

370 

231 


10° 

19 

19 

19 

21 

22 

23 

23 


10 2 

40 

40 

44 

82 

97 

99 

93 


10 4 

42 

42 

47 

94 

220 

241 

181 


10 6 

40 

41 

46 

94 

228 

262 

189 


the lowest order Taylor-Hood element, i.e., (T > 2 ,'Pi,Vi) Lagrange finite elements 
for ( u,ps,p F ). In the experiment, block-diagonal preconditioner B based on the 
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norms in (13.91) . i.e., 

/-A ' 1 0 0 

(3.11) B= 0 I- 1 0 

\ 0 0 (I-nAy 1 

is implemented using FEniCS with Hypre. As above, the exact inverses are replaced 
by suitable algebraic multigrid operators. A heuristic way to validate this choice 
of preconditioner is to consider a special case of (13.41) with A = 1. Noting that 
this special case has the systems in Example 12.11 with A = 1 and in Example 12.21 as 
subsystems, the operator B in (13.111) can be viewed as a combination of the robust 
preconditioners in those examples. 

Numbers of iterations for different constant parameter values are given in Ta¬ 
ble [4] When A is not too large the numbers of iterations are more or less robust 
with respect to variation of n. However, the number of iterations clearly increases 
with increasing A, and the increment is quite large when k is small. The growth 
in iterations as A increases is milder when k is larger, which supports our heuristic 
analysis of “partial remedy”, that (kA7pf,^Pf) may play a role of a stabilization 
term for the L 2 norm of pf ■ 

3.2. A new three-field formulation. The discussion in the previous subsection 
suggests that we need a different formulation to obtain parameter-robust precon¬ 
ditioners. We will present such a new formulation of the system (3T]) here, where 
parameters A, a , k are assumed to be as specified in the beginning of this section. 
In particular, A and k are allowed to be functions of the spatial domain. 

The main obstacle in the discussion above was that the inf-sup condition (13.101) 
is not fulfilled. To circumvent this, we use a different system with unknowns 
(■ u,pt,Pf ), where pr '■= —Adivu + apF- With this new unknown p-r, which will 
be called total pressure, we can rewrite the equations (13.11) as 

- dive (it) + Vp T = /, 

— divit — A -1 (pt ^ ®-Pf) = 0 , 

\~ 1 {apT — ZoPpp) + div(KVp F ) = g. 

The matrix form of this system is 

(u\ /- div e V 0 \ ( u \ ( f \ 

(3.12) A \pr I := —div —A -1 aA -1 PT = 0 . 

\PfJ \ 0 aA -1 — 2 a 2 A _1 + div(«V)/ \PfJ \~9j 

With function spaces Hq, L 2 , Hq for u , pt, pf, respectively, corresponding to 
Dirichlet boundary conditions for u and p F , we obtain a variational form 

(e(n),e(v)) - (divu,p r ) = v e H\, 

(3.13) —(divu,g T ) - (A _ Vt,(?t) + (aA _ 1 p F , <?t) = 0 , q T € L 2 , 

(i a\~ 1 p T ,qF) - 2(ct 2 A~ 1 pF, <7 f) - («VpF, V<7 f) = ( g,qF), qF G H], 

Recall that we used the decomposition p = p m + po and the stabilization term 
in order to obtain the stability of the system in Example 12.21 We need a similar 
argument to show the stability of (13.131) due to the same reason, divi?Q C L 2 . 
Denoting the mean-value zero part of qr by qr, o as in (12.51) . we define norms by 

(3.14) ||e(tt)|| 0 , ((A _ 1 pt,Pt) + Ibr.ollo)' » ((a 2 A _ 1 pF ,Pf) + (kVp f , Vp F )) 5 
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for u £ Hi, p T € L 2 , pf £ H g. Let us denote these spaces with the norms in 
(Id. 1 dll by V , Qt, Qf, and let X = V x Qt x Q f . Then it can be shown that all 
the bilinear forms in (13.131) are uniformly bounded in A, a, k. In other words, the 
operator A appearing in (13.1211 is a bounded linear map from X to X* and its norm 
is independent of the three parameters. Here the norm on the space X* is derived 
from the norm on X exactly as we explained in Example 12.11 

The following theorem implies that A is invertible and that the inverse is a map 
from X* to X with operator norm independent of the three parameters. 

Theorem 3.2. For the system (13.131) there exists a constant p > 0, independent 
of X, a, k satisfying (13.21) . such that the following inf-sup condition holds: 

. , (A(u,pt,Pf), {v,q T ,q F ))( X . X ) 

ln f Slip — - -r—r- - — - > P- 

(u,PT,pf)£X q T ,q F )ex \\{u,pt,Pf)\\x\\(v, qr, qF)\\x 

Proof. To prove the inf-sup condition, we will use a standard technique: For given 
(0,0,0) ^ (u,pt,Pf) £ X, we will find ( v,qT,qF ) £ X such that 

(3-15) \\(v,q T ,qF)\\x < Ci\\(u,p t ,Pf)\\x, 

(3-16) (A{u,p T ,p F ),(v,q T ,q F ))(x* ,x) > C 2 \\(u,Pt,Pf)\\%i 

with positive constants C\, C 2 which are independent of A, re, and a. From these 
two inequalities we obtain that the desired inf-sup condition holds with p = C 2 /C 1 . 

Suppose that (0,0,0) ^ (u,pt,Pf) £ X is given. Recall that p F ,o is the mean- 
value zero part of pt ■ It is well-known from the theory of Stokes equation, cf. |39] 
Theorem 5.1] that there exists a constant Po > 0, depending only on the domain 
12, and w £ V, such that 

(3.17) (div w,pt) = |br,o||o> (e(w), e(w)) < ^oIIpt.oIIo- 

We set v = u — Sow, q F = —Pt, qF = —pf with a constant Sq which will be 
determined later. One can check that 

\\{v,qT,q F )\\x < \j2(l + 5lPl)\\{u,pT,p F )\\ x , 

and (13.151) follows, if is independent of the parameters of our interest. 

To establish (13.161) and determine So, we use the chosen v , q F , qF, and (13.171) 
to obtain 


(• A{u,p T ,p F ), (v, q T , q F ))(x*,x) 

= ll«llv-^o(e(M),e(^)) + <5 0 (divm,p T ) 

(3.18) + ((X~ 1 Pt,Pt) + 2(a 2 \~ 1 p F ,p F ) - 2(a\~ 1 p T ,p F )) + (reVp F , X7p F ) 

= ||«||v- 5 o(e(«), e ( w )) + <5o||pt,o||o 
+ ((A ~ 1 p Tl pr) + 2(a 2 X~ 1 p F ,p F ) - 2(a\~ 1 p T ,p F )) + (reVp F , Vp F )- 
By Young’s inequality and (13.171) , we also have 


{e(u),e{w)) < + y IMIv < y^lMlv + ^y 

Using the above inequality with the choice do = <5o = Pq 2 , 
(3.19) |M|^-£ 0 (e(u), £(«;)) + ^o||pt,o||o > \\H\v X 


-\\pt,o\\ 2 o, 


we derive 

y Ibr.ollo- 


V6» 0 > 0. 
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Again by Young’s inequality, for any > 0, 

2(aX~ 1 p T ,PF) = 2(\~ 1/2 p T ,a\~ 1/2 p F ) < 2d 1 (X ~ 1 p T ,pr) + -^-(a 2 X~ 1 pf,Pf)- 

Zu\ 

If we take = 3/8, then we get 

(3.20) {(X~ 1 Pt,Pt) + 2(a 2 X 1 pf,Pf) — 2(aA~ 1 p T ,PF)) 

1 2 

- ^(^"Vt^t) + ^(« 2 A~Vf,Pf). 

The inequality (13.161) is obtained by combining (13.181) . (|3.19l) . and (13.201) . Finally, 
we remark that the constants C\ and Ci in (13.151) (13.161) depend only on Sq, which 
is /3q 2 in the argument, so they are independent of A, a, n. □ 


Remark 3.3. The set up in Theorem 13.21 is suitable for homogeneous Dirichlet 
boundary conditions for the displacement u and the fluid pressure p F - However, 
similar result can be obtained for more general boundary conditions. For this, we 
first review possible boundary conditions for Biot’s model. Suppose that there are 
two partitions of dFl, 

(3.21) dn=r p uT f , dn = r d uv t7 

with |r p |, |r d | > 0, i.e., the Lebesgue measures of T p and r d are positive. General 
homogeneous boundary conditions of m are given by 

p F (t) = 0 on r p , —KWp F {t) • n = 0 onT/, 
u(t) =0 on r d , S.(t)n = 0 on r t , 

for time variable t G [0,T], T > 0, in which n is the outward unit normal vector 
field on dFl and <r(t ) := 2 pe[u(t)) + (Adivit(f) — ap^(t))/. The conditions for 
Pf is a combination of pressure-flux boundary condition as in Darcy flow and the 
conditions for it is a combination of displacement-traction boundary conditions as 
in elasticity problems. The proper function spaces for this variational formulation 
are 


(3.22) y := {« £ H 1 : u|r d = 0}, Qt = L 2 , Qf := {qF £ H 1 : gF|r p = 0}, 
for u, PT-, and p F - When ^ dfl, we choose parameter-dependent norms by 

(3.23) ||e(u)||o, ||pt||o, {{a 2 X~ l p Fl p F ) + (kVpf, Vp F )) 2 ■ 


We can prove a stability result similar to Theorem 13.21 with these norms. In fact, 
the proof is easier in this case since the inf-sup condition 


(3.24) 


. (divu,p T ) 

mi , sup F-rrj, -rp > do 

PTSL 2 „ ef fi ||v||i||pt||o 


holds, and therefore a decomposition of pr, into its mean-value and mean-value 
zero components, is not necessary. We omit the details since the proof is completely 
analogous to the proof of Theorem 13.21 
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4. Discretization and construction of preconditioners 


In this section we propose finite element discretizations of the three field formu¬ 
lation introduced above, and show that parameter-robust preconditioners for the 
discrete problems can be found. In contrast to the discussion above, we will first 
consider problems with general boundary conditions, i.e., boundary conditions with 
Td ^ dCl, cf. (111.211) . The reason for this reversed order is that the construction 
of preconditioners in the case when T^ = dQ, resulting in the choice V = Hq, 
requires a nontrivial technical discussion. 

We have shown above that (jli.llll) is a linear system with parameter-robust 
stability for the function spaces V, Qt , Qf with parameter-dependent norms 
given by (13.141) . If we discretize the system (13.131) with the finite element spaces 
V h C , QT,h tZ bf t , LfF,h td f , then the discrete counterpart of (|fi_4_il|) is to find 
( Uh:PT,h:PF,h ) £Vh,x Qr,h x QF,h such that 

(e(« h),§.(v)) - (di vv,p T ,h) = Vv £ V h , 

(4.1) — (divii/i, qr) — (A 1 pr,h, <?t) + (<*A 1 pF,h,qr) = 0, Vqr £ Qr,h, 

(aA ~ 1 pT,h, qF ) — 2(o; 2 A~Vf,?i, q F ) — {XS7pF,h,Vqp) = ( g , <?f), £ Qf,7i- 


A basic stability assumption for this discretization is that the pair Vh x Qr,h 
satisfies a discrete version of (13.241) . i.e., 


(4.2) 


inf sup 


(divt>,p T ) 


PT eQ T . hv 7v h IMIillHIo 


> A) > 0, 


where fio is independent of h. In other words, Vh x Qr,h is a stable Stokes pair. 


Theorem 4.1. Suppose that V, Qt, Qf ore as in (13.221) with T^ ^ <9f l and T p 
as in (13.211) , and that Vh C V, Qtm C Qt, QF,h C Qf ore corresponding finite 
element spaces. Furthermore, assume that the pair Vh x Qr,h satisfies the inf-sup 
condition m- Let Xh = Vh x Qtm x QF,h be the Hilbert space with norm given in 
(13.231) . and Ah ■ Xh —> X£ the operator given by (|4.1I) . Then there exists a constant 
/? > 0, such that 


(4.3) 


inf sup 

(u,p T ,p F )GX h (v,q T ,q F )eX h 


( Ah(u,p T ,PF ), ( v ,qT,qF))(X*,X h ) 
W(u,PT,PF)\\x h \\{v,q T ,qF)\\x h 


>d: 


for all parameters A, a, and n satisfying 


We do not prove this result here since the proof is completely analogous to the 
proof of Theorem 13.21 We observe that the norms given in (13.141) shows that a 
preconditioner of the form 


/-A ” 1 0 0 

(4.4) B= 0 I~ l 0 

V 0 0 (a 2 A -1 / — div(ttV )) _1 

will be a parameter-robust preconditioner. 

We now turn to the case with T^ = dLl such that V = H^. We recall that 
Lq is the space of L 2 functions with mean value zero. The proper discrete inf-sup 
condition in this case takes the form 


inf sup 

PT&QT.hFLl vGV h 


(div v,p T ) 


lIIt 


T || 0 


> A, > o, 


(4.5) 
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where again f3o is independent of h. 

The following is a discrete analogue of Theorem l3.2l and its proof is completely 
analogous to the proof of that theorem. 

Theorem 4.2. Suppose that V, Qt, Qf are as in (13.221) with T^ = dCl and T p 
as in ([3.211) . and that V h C V, Qt,h C Qt, Q F,h C Qf are corresponding finite 
element spaces. Furthermore, assume that the pair V/, x Qr,h satisfy the inf-sup 
condition (fop . Let Xh = Vh x Qr,h x Qf.h be the Hilbert space with norm given 
in (|3.14l) . and Ah '■ Xh —> Xj* the operator given by (14.11) . Then there is a constant 
ft > 0 such that 631 ) holds for all parameters A, a, and k satisfying (] 3.2 1) . 


There exist a number of choices of stable Stokes pairs Vh x Qr,h , and in Section 
6 below we will present numerical results for two examples, the lowest order Taylor- 
Hood element and the MINI element. For more examples of stable Stokes pairs we 
refer to [29] [39]. The parameter-dependent norms in (13.141) suggest a block diagonal 
preconditioner of the form 


(4.6) 


B = 


t- A- 
0 
0 


0 

(A^J + Io) 
0 


-l 


(a 2 A 1 I — div(«V)) 


for the continuous system. We recall that / is the Riesz map of Qt into its dual 
Qt, and Iq is the corresponding map into the dual of Qt D Lq. The first and third 
blocks of this block diagonal operator are inverses of standard second-order elliptic 
operators, and corresponding preconditioners to replace the exact inverses in the 
discrete case are well-studied. In contrast, the operator in the second block is less 
standard, and, as far as we know, a construction of an effective preconditioner to 
replace it has not been proposed. We will discuss such a construction below. 


5. A PRECONDITIONER FOR THE OPERATOR A + Iq 

Throughout this section the parameter A is assumed to be a constant. We recall 
from the discussion above that in order to construct an effective block diagonal 
preconditioner of the form (14.61) we need to replace the inverse of the operator 
A" 1 / + Iq by a spectrally equivalent operator which can be cheaply evaluated. In 
fact, when A > 1 the operators A -1 / + Iq and A -1 / m + Iq are spectrally equivalent, 
so it is enough to approximate the inverse of the latter. 

Let N be the dimension of Qr,h and i be the standard nodal basis of 

QT,h- Let Iq be the constant function on Q with value 1/^/fsTf where |f2| is the 
volume of Cl. Denoting the mean-value zero and mean-value parts of <f l by <f >q and 
(j) l m as before, we can observe that 

(5.1) 4>m = miln, to* := (<(>\ Iq), = </>(, + VI <i<N. 

If we let M, Mo, M m be mass matrices corresponding to the operators J, Iq, and 
I m , then their (i, j)-entries are 

M 0 »M*, j) := (<fi, fa), VI <i,j<N, 

and the matrix corresponding to A ~ 1 I m + A) is A -1 M m + M 0 . Since (<,b l ,(jS) = 
(4 >q + 4> l m , <f>l + (fir ,) = ((fo , ) + {<f> l m , (fin) , one can see M = M 0 + M m , and therefore 
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A 1 M m + Mo = M + (A 1 — l)M m . In addition, observe that M m (i, j) = rmmj by 
(15.11) , so 

(5.2) A _ 1 M m +M 0 =M+ (A -1 - l)mm T , 


with 


(5.3) 


m = 


(m 1 \ 

m 2 


\m N J 


To construct a preconditioner we need to find an approximate inverse of A 'M m + 
Mo- Since M is positive definite, we can rewrite the right-hand side of (15.21) as 

(5.4) Ma := M + (A ” 1 — l)mm T = (I + (A -1 — l)mm T M _ 1 )M, 


where I is the N x N identity matrix. Recall the Sherman-Morrison-Woodbury 
formula, 


T 

(5.5) (I + uv T ) _1 = I - =— , u, v £ M. n with u T v ^ — 1. 

v ' v ’ l + u T v 

We will use it to find the inverse of I + (c -1 — l)mm T i _1 for a constant c ^ 0. 

Lemma 5 . 1 . Let w = (1 • • • 1) T £ R^. For the mass matrix M and m in (15.31) . 
the following two identities holds: 

(5.6) Mw = i/|n|m, m T w = \/|f2|. 


Proof. Note that 
(5 ' 7> 

because {<jF}i.<j<N is the standard nodal basis and no boundary condition is im¬ 
posed on Qr,h- 

If we consider the i-th row of the left-hand side of the first identity in (15.61) . then 
the definition of M, & and m give 

N N / N \ 

E M (*,i) = Ew^') = ^>E^ = VW\mi. 

3= 1 3 = 1 V i =1 / 

This proves the first identity in (15.61) . The second identity follows by 

N N 

5> = E(^> 1q ) = Vl^Rln, In) = VW\- 

i=1 i=1 

□ 


Corollary 5.2. For the mass matrix M, m in (|5.3|1 . w in Lemma 1,5.11 and any 
constant c / 0, the following holds: 

(5.8) (I + (c -1 — l)mm T M _ y l = I + (c — l)(-\/i^) - 1 mvv T . 

Proof. Since M is symmetric, the first identity in (15.61) gives m T M" - 1 = (,/jnj )- 1 W T . 
If we set u = (c _1 — l)m and v T = m r M _1 , then the second identity in (15.61) gives 
u T v = c -1 — 1^—1. The assertion follows from (15.51) . □ 
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Theorem 5.3. Let Y\ = (I + amw r ) 1 with a = (—1 + y/X)/y/\Q\. Then, for 
Ma in m , Ma = VaMVJ. Thus, ?/B is a preconditioner of M with condition 
number K, thenY^WVf 1 is a preconditioner of Ma with same condition number. 

Proof. From the definition of Va and the second identity in (15.61) . we can see 

V^ 2 = (I + amw T ) 2 = I + (2a + a 2 \/jnj')mw T = I + (A — l)(-y/jn[) _1 mw T . 

By the identity (15.81) with c = A and (15.41) , we have Ma = V^M. If we use (15.81) for 
V^ 1 , then one can verify that Va = I + amm T M -1 with a = 1 — (a/A) - 1 . From 
this expression of Va, it is easy to check that VaM = MVj, so Ma = V|M = 
VaMV^. The assertion for preconditioner Vi^BV^ 1 follows from the identity 
MaV" t BV" 1 M a = VaMBMVJ. □ 

For the preconditioner D for M, it is known that the Jacobi preconditioner, i.e., 
the inverse of diagonal of mass matrix as a preconditioner has explicit condition 
number bounds |40| . If Qr,h is the piecewise linear continuous finite element, then 
the Jacobi preconditioner D for M is a constant multiple of the diagonal matrix 
diag(m]" 1 , mif 1 ,..., mf^). As a consequence, wm T B = Omw T , so V^BV^ 1 can 
be reduced to UV^ 2 = B(I + (A — l)(A/|Ti[)" 1 mw T ). 

There is one caution in the implementation of the preconditioner Vi^BV^ 1 
because mw T and mm 2 in V^ 1 and Ma are dense matrices in general. We re¬ 
mark that the minimum residual method requires only matrix-vector multiplica¬ 
tion operations. Therefore, to avoid computation with these dense matrices, we 
use the structure of the matrix mw T . More precisely, mw T v for an R w -vector 
v can be computed with two operations, the inner product w T v and constant- 
vector multiplication (w T v)m. Similarly, we can avoid generating mm T in Ma = 
M + (A -1 — l)mm T . Finally, we remark that the preconditioner Vi^BV^ 1 is useful 
when a piecewise discontinuous finite element is used for Qr,h because mw T and 
Ma are not sparse. 

Table 5. Boundary conditions (BC), preconditioners (PC), and 
finite elements of test cases. The first three cases use the lowest or¬ 
der Taylor- Hood element and the last case uses the MINI element 
(B = vector-valued bubble function). 




BC 


PC 

finite elements 

system size for N 

32 64 128 

Case 1 

r d y dn, r p 

= an 

(l4~4l> 

V2-V1-V1 

10628 

41732 

165380 

Case 2 

r d 

= on, r p 

= an 

(l4~6ll 

V2-V1-V1 

10628 

41732 

165380 

Case 3 

r d 

= an, r p 

= an 

m 

'P2-P1-V1 

10628 

41732 

165380 

Case 4 

r d + on, r P 

= an 

m 

(Vi + B)-Vi-Vi 

8452 

33284 

132100 


6. Numerical results 

In this section we present some numerical results which illustrate our theoretical 
results for the proposed preconditioners ( 14 . 41 ) and ( 14 . 61 ) . As before, all numerical 
experiments are carried out using FEniCS with Hypre algebraic multigrid operators 
as replacements for the exact inverses appearing in the first and third block of ( 14 . 41 ) 
and ( 14 . 61 ) . Furthermore, in the preconditioner of the form ( 14 . 61 ) , the second block is 
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TABLE 6. Numbers of iteration and condition numbers of Case 1 (cf. 
Table[5]). (fl = unit square, partitioned as bisections of NxN rectangles, 
convergence criterion with relative residual of 10 -6 ) 


N 

a 

A 

10° 

IQ - 4 

10“ 8 

10“ 12 



10° 

33 (3.8) 

43 (6.3) 

47 (7.6) 

47 (7.6) 


10° 

10 4 

52 (21.7) 

52 (21.7) 

65 (21.7) 

63 (21.7) 



10 8 

52 (21.7) 

54 (21.7) 

52 (21.7) 

62 (21.7) 



10° 

33 (3.8) 

33 (3.9) 

43 (6.3) 

47 (7.6) 

32 

1(T 2 

10 4 

52 (21.7) 

52 (21.7) 

52 (21.7) 

63 (21.7) 



10 8 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 



10 u 

33 (3.8) 

33 (3.8) 

33 (3.8) 

43 (6.3) 


10" 4 

10 4 

50 (21.7) 

52 (21.7) 

50 (21.7) 

52 (21.7) 



10 8 

54 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 



10° 

33 (3.9) 

40 (5.6) 

47 (7.6) 

47 (7.6) 


10° 

10 4 

52 (21.7) 

52 (21.7) 

63 (21.7) 

63 (21.7) 



10 8 

46 (21.7) 

52 (21.7) 

52 (21.7) 

62 (21.7) 



10 u 

33 (3.8) 

33 (3.8) 

40 (5.6) 

47 (7.6) 

64 

1(T 2 

10 4 

52 (21.7) 

52 (21.7) 

52 (21.7) 

58 (21.7) 



10 s 

52 (21.7) 

46 (21.7) 

52 (21.7) 

48 (21.7) 



10 u 

33 (3.9) 

33 (3.8) 

33 (3.8) 

40 (5.6) 


10" 4 

10 4 

50 (21.7) 

46 (21.7) 

50 (21.7) 

52 (21.7) 



10 8 

52 (21.7) 

52 (21.7) 

50 (21.7) 

52 (21.7) 



10° 

32 (3.8) 

39 (5.4) 

46 (7.7) 

45 (7.7) 


10° 

10 4 

51 (21.7) 

52 (21.7) 

61 (21.7) 

59 (21.7) 



10 8 

51 (21.7) 

50 (21.7) 

50 (21.7) 

54 (21.7) 



10 u 

33 (3.8) 

33 (3.8) 

39 (5.2) 

46 (7.7) 

128 

10 -2 

10 4 

48 (21.7) 

45 (21.7) 

52 (21.7) 

58 (21.7) 



10 8 

50 (21.7) 

52 (21.7) 

50 (21.7) 

52 (21.7) 



10° 

33 (3.8) 

32 (3.8) 

33 (3.8) 

39 (5.2) 


10" 4 

10 4 

44 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 



10 8 

50 (21.7) 

48 (21.7) 

48 (21.7) 

50 (21.7) 


constructed by using the technique outlined in Section O while the standard Jacobi 
preconditioner is used in the second block of (14.411 . 

In all the experiments the domain fl is unit square. We show numerical results 
for four different combinations of boundary conditions, finite element spaces, and 
preconditioners. The different combinations are presented as Case 1-4 in Table [5] 
In Case 1 and 4 the statement Tj ^ dtt means that is taken as in Example 
3.1, while problems with = <9f2 are consider in Case 2 and Case 3. We compare 
numerical results obtained by the two preconditioners with structure of the form 
gUl) and CEO). The result of Theorem ro suggests that preconditioners of the 
form (14.61) are more robust than the ones of the form (14.41) in the case of Dirichlet 
boundary conditions. In other words, we expect that the results of Case 2 are 
more robust than the ones of Case 3. However, the system preconditioned with 
a preconditioner of the form (14.41) has only one bad eigenvalue. Therefore, as in 
Example 12.21 we can expect small differences in the number of iterations. Finally, 


















Table 7. Numbers of iteration of test cases in Table[S] (f2 = unit 
square, partitioned into bisections of N xN rectangles, convergence 
criterion with relative residual of 10~ 6 ) 
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Case 1 


Case 2 


Case 3 


Case 4 





N 



N 



N 



N 


K, 

a 

A 

32 

64 

128 

32 

64 

128 

32 

64 

128 

32 

64 

128 



10° 

33 

33 

32 

29 

29 

29 

29 

29 

29 

34 

34 

34 


10° 

10 4 

52 

52 

51 

46 

46 

46 

66 

45 

44 

60 

61 

60 

10° 


10 8 

52 

46 

51 

46 

46 

45 

44 

44 

43 

60 

61 

60 


10 u 

33 

33 

33 

29 

29 

29 

29 

29 

29 

34 

34 

34 


10" 4 

10 4 

50 

50 

44 

46 

46 

44 

66 

44 

44 

60 

60 

60 



10 s 

54 

52 

50 

46 

46 

45 

45 

44 

43 

60 

60 

60 



10 u 

43 

40 

39 

39 

38 

36 

39 

38 

36 

46 

43 

40 


10° 

10 4 

52 

52 

52 

46 

46 

45 

68 

44 

44 

60 

61 

60 

1(T 4 


10 8 

54 

52 

50 

46 

46 

45 

66 

44 

43 

60 

60 

60 


10 u 

33 

33 

32 

29 

29 

29 

29 

29 

29 

34 

34 

34 


10“ 4 

10 4 

52 

46 

52 

46 

46 

45 

66 

44 

42 

60 

60 

60 



10 s 

52 

52 

48 

46 

46 

45 

44 

44 

43 

60 

61 

60 



10° 

47 

47 

46 

42 

42 

42 

42 

42 

42 

52 

52 

52 


10° 

10 4 

65 

63 

61 

61 

59 

58 

60 

58 

57 

73 

73 

72 

10~ 8 


10 8 

52 

52 

50 

46 

46 

45 

44 

44 

43 

60 

60 

60 


10 u 

47 

33 

33 

29 

29 

29 

29 

29 

29 

34 

34 

34 


10“ 4 

10 4 

65 

50 

52 

46 

45 

45 

66 

44 

43 

60 

60 

60 



10 8 

52 

50 

48 

46 

46 

44 

44 

44 

40 

60 

60 

60 



10 u 

47 

47 

45 

42 

42 

42 

42 

42 

42 

52 

52 

52 


10° 

10 4 

63 

63 

59 

58 

58 

57 

57 

56 

56 

72 

72 

72 

10 -12 


10 s 

62 

62 

54 

58 

58 

50 

57 

56 

54 

72 

72 

71 


10 u 

43 

40 

39 

39 

38 

36 

37 

38 

36 

47 

43 

40 


10“ 4 

10 4 

52 

52 

52 

46 

46 

44 

66 

44 

44 

60 

61 

60 



10 8 

52 

52 

50 

46 

45 

44 

44 

43 

43 

60 

61 

60 


in Case 4 we use the MINI element instead of the Taylor-Hood element in order to 
show that our results are robust with respect to the choice of finite element spaces, 
as long as they fulfill the assumptions of the theory. In most of the examples the 
parameters A, a and k are taken to be constants. However, in the last experiment, 
presented in Tabic [HI k varies with the spatial variable. 

In Table [6l we present numbers of iteration of Case 1. The results are fairly 
robust with respect to parameter changes and mesh refinements. To compare ro¬ 
bustness of preconditioners of all the cases, we present numbers of iteration for all 
the different four cases in Table [7] As expected, the results of Case 2 are slightly 
better than the ones of Case 3, in particular for N = 32. Although there are no 
remarkable differences in the presented results, in the full numerical results which 
are not included here, the results for Case 3 shows that this method sometimes need 
about 30 — 45% more iterations than those of Case 2. In Case 4 we use the MINI 
element instead of the Taylor-Hood element. Although the numbers of iteration 
are larger than those for the Taylor-Hood element, the results are still quite robust 
with respect to changes of parameters. As the final experiment, a model problem 
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TABLE 8. Numbers of iteration and condition numbers of Case 1 (cf. 
Table [5]) with nonconstant k. (fl = unit square, partitioned as bisections 
of N x N rectangles, fli = {(x,y) : 0 < x < 1,1/4 < y < 3/4}, k = 1 
on Q \ fli, convergence criterion with relative residual of 10 -6 ) 







k on 



N 

a 

A 

10” 2 

10" 4 

10“ 6 

10" 8 

10- 10 



10° 

34 (4.5) 

42 (6.1) 

46 (7.4) 

46 (7.5) 

46 (7.4) 


10° 

10 4 

54 (21.7) 

52 (21.7) 

56 (21.7) 

65 (21.7) 

61 (21.7) 



10 8 

54 (21.7) 

54 (21.7) 

46 (21.7) 

52 (21.7) 

55 (21.7) 



10 u 

33 (3.8) 

33 (3.8) 

34 (4.4) 

42 (6.1) 

46 (7.4) 

32 

10 -2 

10 4 

52 (21.7) 

48 (21.7) 

50 (21.7) 

53 (21.7) 

56 (21.7) 



10 8 

52 (21.7) 

52 (21.7) 

50 (21.7) 

52 (21.7) 

46 (21.7) 



10° 

33 (3.8) 

33 (3.8) 

33 (3.8) 

33 (3.8) 

34 (4.5) 


10" 4 

10 4 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 



10 8 

52 (21.7) 

52 (21.7) 

50 (21.7) 

53 (21.7) 

50 (21.7) 



10 u 

34 (4.5) 

39 (5.4) 

44 (7.4) 

46 (7.5) 

46 (7.5) 


10° 

10 4 

52 (21.7) 

52 (21.7) 

55 (21.7) 

63 (21.7) 

60 (21.7) 



10 8 

52 (21.7) 

46 (21.7) 

46 (21.7) 

50 (21.7) 

51 (21.7) 



10° 

33 (3.8) 

33 (3.8) 

34 (4.4) 

39 (5.4) 

45 (7.4) 

64 

10“ 2 

10 4 

46 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 

51 (21.7) 



10 8 

48 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 



10 u 

33 (3.8) 

33 (3.8) 

33 (3.8) 

33 (3.8) 

34 (4.4) 


10" 4 

10 4 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 



10 8 

52 (21.7) 

52 (21.7) 

52 (21.7) 

52 (21.7) 

46 (21.7) 



10 u 

34 (4.4) 

37 (5.0) 

44 (7.2) 

46 (7.5) 

44 (7.5) 


10° 

10 4 

50 (21.7) 

48 (21.7) 

49 (21.7) 

56 (21.7) 

61 (21.7) 



10 8 

50 (21.7) 

50 (21.7) 

52 (21.7) 

50 (21.7) 

47 (21.7) 



10 u 

33 (3.8) 

32 (3.9) 

34 (4.4) 

37 (5.0) 

44 (7.1) 

128 

10 -2 

10 4 

50 (21.7) 

46 (21.7) 

50 (21.7) 

52 (21.7) 

49 (21.7) 



10 s 

52 (21.7) 

52 (21.7) 

50 (21.7) 

52 (21.7) 

50 (21.7) 



10 u 

32 (3.8) 

33 (3.9) 

33 (3.8) 

33 (3.9) 

34 (4.4) 


10" 4 

10 4 

52 (21.7) 

50 (21.7) 

52 (21.7) 

49 (21.7) 

46 (21.7) 



10 8 

52 (21.7) 

50 (21.7) 

48 (21.7) 

52 (21.7) 

50 (21.7) 


with nonconstant k is considered. We assume that k is small on 

Sdi = {(x, y ) : 0 < x < 1,1/4 < y < 3/4} C ff, 

and k = 1 on f] \ The numerical results in Table [5] are fairly robust for mesh 
refinements and changes of parameters, including high contrasts of k. 

7. Conclusion 

We have studied parameter-robust discretizations and construction of precon¬ 
ditioners for Biot’s consolidation model. To apply the framework of |25j we have 
proposed a new three-field formulation of the Biot system. We have showed that 
preconditioners based on mapping properties and parameter-dependent norms are 
robust with respect to variations of the model parameters, choice of finite element 
spaces satisfying the proper stability condition, and the discretization parameters. 
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In particular, the variations of parameters in our consideration cover large shear 
and bulk elastic moduli, small hydraulic conductivity, small time-step, including 
the ranges of interest in geophysics and computational biomechanics applications. 
Furthermore, our theoretical results are confirmed by a number of numerical ex¬ 
periments. 


References 

[1] R. E. Showalter. Diffusion in poro-elastic media. J. Math. Anal. Appl., 251(l):310-340, 2000. 

[2] M. B. Reed. An investigation of numerical errors in the analysis of consolidation by finite 
elements. Internat. J. Numer. Analyt. Methods Geomech., 8(3):243-257, 1984. 

[3] P. A. Vermeer and A. Verruijt. An accuracy condition for consolidation by finite elements. 
Internat. J. Numer. Analyt. Methods Geomech., 5(1):1-14, 1981. 

[4] O. C. Zienkiewicz and T. Shiomi. Dynamic behaviour of saturated porous media; the gen¬ 
eralized Biot formulation and its numerical solution. Internat. J. Numer. Analyt. Methods 
Geomech., 8(1):71—96, 1984. 

[5] M. A. Murad and A. F. D. Loula. Improved accuracy in finite element analysis of Biot’s 
consolidation problem. Comput. Methods Appl. Mech. Engrg., 95(3):359-382, 1992. 

[6] M. A. Murad and A. F. D. Loula. On stability and convergence of finite element approxima¬ 
tions of Biot’s consolidation problem. Internat. J. Numer. Methods Engrg., 37(4):645-667, 
1994. 

[7] M. A. Murad, V. Thomee, and A. F. D. Loula. Asymptotic behavior of semidiscrete finite- 
element approximations of Biot’s consolidation problem. SIAM J. Numer. Anal., 33(3): 1065- 
1083, 1996. 

[8] J. Korsawe and G. Starke. A least-squares mixed finite element method for Biot’s consolida¬ 
tion problem in porous media. SIAM J. Numer. Anal., 43(1):318-339, 2005. 

[9] Y. Chen, Y. Luo, and M. Feng. Analysis of a discontinuous Galerkin method for the Biot’s 
consolidation problem. Appl. Math. Comput., 219(17):9043-9056, 2013. 

[10] P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin finite element 
methods for poroelasticity. I. The continuous in time case. Comput. Geosci., 11 (2): 131-144, 
2007. 

[11] P. J. Phillips and M. F. Wheeler. A coupling of mixed and continuous Galerkin finite element 
methods for poroelasticity. II. The discrete-in-time case. Comput. Geosci., 11(2):145-158, 
2007. 

[12] P. J. Phillips and M. F. Wheeler. A coupling of mixed and discontinuous Galerkin finite- 
element methods for poroelasticity. Comput. Geosci., 12(4):417-435, 2008. 

[13] S.-Y. Yi. A coupling of nonconforming and mixed finite element methods for Biot’s consoli¬ 
dation model. Numer. Methods Partial Differ. Equ., 29(5): 1749-1777, 2013. 

[14] F. Ben-Hatira, K. Saidane, and A. Mrabet. A finite element modeling of the human lumbar 
unit including the spinal cord. J. Biomed. Sci. Eng., 5:146-152, 2012. 

[15] J. H. Smith and J. A. Humphrey. Interstitial transport and transvascular fluid exchange 
during infusion into brain and tumor tissue. Microvasc. Res., 73(l):58-73, 2007. 

[16] K. Stpverud, M. Alnaes, H. P. Langtangen, V. Haughton, and K.-A. Mardal. Poroelastic 
modeling of syringomyelia - a systematic study of the effects of pia mater, central canal, 
median fissure, white and grey matter on pressure wave propagation and fluid movement 
within the cervical spinal cord. Comput. Methods Biomech. Biomed. Engin. Accepted 2015. 

[17] O. Coussy. Poromechanics. John Wiley & Sons, 2004. 

[18] H. F. Wang. Theory of linear poroelasticity. Princeton Series in Geophysics, Princeton Uni¬ 
versity Press, Princeton, NJ, 2000. 

[19] O. Axelsson, R. Blaheta, and P. Byczanski. Stable discretization of poroelasticity prob¬ 
lems and efficient preconditioners for arising saddle point type matrices. Comput. Vis. Sci., 
15(4):191-207, 2012. 

[20] S.-H. Chan, K.-K. Phoon, and F. H. Lee. A modified Jacobi preconditioner for solving ill- 
conditioned Biot’s consolidation equations using symmetric quasi-minimal residual method. 
Internat. J. Numer. Analyt. Methods Geomech., 25(10): 1001-1025, 2001. 


24 


JEONGHUN J. LEE, KENT-ANDRE MARDAL AND RAGNAR WINTHER 


[21] J. B. Haga, H. Osnes, and H. P. Langtangen. A parallel block preconditioner for large-scale 
poroelasticity with highly heterogeneous material parameters. Comput. Geosci., 16(3):723- 
734, 2012. 

[22] K. K. Phoon, K. C. Toh, S. H. Chan, and F. H. Lee. An efficient diagonal preconditioner 
for finite element solution of Biot’s consolidation equations. Int. J. Numer. Methods Eng., 
55(4):377-400, 2002. 

[23] S. Rhebergen, G. N. Wells, R. F. Katz, and A. J. Wathen. Analysis of block preconditioners 
for models of coupled magma/mantle dynamics. SIAM J. Sci. Comput ., 36(4):A1960-A1977, 

2014. 

[24] S. Rhebergen, G. N. Wells, A. J. Wathen, and R. F. Katz. Three-field block preconditioners 
for models of coupled magma/mantle dynamics. SIAM J. Sci. Comput ., 37(5):A2270-A2294, 

2015. 

[25] K.-A. Mardal and R. Winther. Preconditioning discretizations of systems of partial differential 
equations. Numer. Linear Algebra Appl., 18(1): 1-40, 2011. 

[26] R. D. Falgout and U. M. Yang. Computational Science -— ICCS 2002: International Confer¬ 
ence Amsterdam, The Netherlands, April 21-24, 2002 Proceedings, Part III, chapter hypre: 
A Library of High Performance Preconditioners, pages 632-641. Springer Berlin Heidelberg, 
Berlin, Heidelberg, 2002. 

[27] K.-A. Mardal and J. B. Haga. Block preconditioning of systems of PDEs. In Automated 
Solution of Differential Equations by the Finite Element Method, pages 643-655. Springer 
Berlin Heidelberg, 2012. 

[28] A. Logg, K.-A. Mardal, and G. N. Wells, editors. Automated solution of differential equations 
by the finite element method, volume 84 of Lecture Notes in Computational Science and 
Engineering. Springer, Heidelberg, 2012. 

[29] F. Brezzi and M. Fortin. Mixed and Hybrid Finite Element Methods , volume 15 of Springer 
Series in computational Mathematics. Springer, 1992. 

[30] O. Axelsson and G. Lindskog. On the rate of convergence of the preconditioned conjugate 
gradient method. Numer. Math., 48(5):499-523, 1986. 

[31] B. F. Nielsen and K.-A. Mardal. Analysis of the minimal residual method applied to ill posed 
optimality systems. SIAM J. Sci. Comput., 35(2):A785-A814, 2013. 

[32] A. Anandarajah. Computational Methods in Elasticity and Plasticity: Solids and Porous 
Media. Springer New York, 2010. 

[33] I. Babuska and M. Suri. Locking effects in the finite element approximation of elasticity 
problems. Numer. Math., 62(4):439-463, 1992. 

[34] D. S. Malkus and T. J. Hughes. Mixed finite element methods reduced and selective integra¬ 
tion techniques: A unification of concepts. Comput. Methods Appl. Mech. Engrg., 15(1):63 - 
81, 1978. 

[35] D. Boffi, F. Brezzi, and M. Fortin. Finite elements for the Stokes problem. In D. Boffi and 
L. Gastaldi, editors, Mixed Finite Elements: Compatibility Conditions, volume 1939 of Lec¬ 
ture Notes in Mathematics. Springer, 2008. 

[36] S. C. Brenner and L.-Y. Sung. Linear finite element methods for planar linear elasticity. Math. 
Comp., 59(200):321-338, 1992. 

[37] M. Crouzeix and P.-A. Raviart. Conforming and nonconforming finite element methods for 
solving the stationary Stokes equations. I. Rev. Frangaise Automat. Informat. Recherche 
Operationnelle Ser. Rouge, 7(R-3):33-75, 1973. 

[38] K.-A. Mardal, X.-C. Tai, and R. Winther. A robust finite element method for Darcy-Stokes 
flow. SIAM J. Numer. Anal., 40(5): 1605-1631, 2002. 

[39] V. Girault and P.-A. Raviart. Finite element methods for Navier-Stokes equations, volume 5 
of Springer Series in Computational Mathematics. Springer-Verlag, Berlin, 1986. 

[40] A. J. Wathen. Realistic eigenvalue bounds for the Galerkin mass matrix. IMA J. Numer. 
Anal, 7(4):449-457, 1987. 


